﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SystemAnalizer.Algorithms
{
   public class LR
    {
        private double[, ,] Discret_matrix;
        private int system_size;
        private int discrete_num;
        private List<double[, ,]> A_matrix;
        private List<double[, ,]> L_matrix;
        private List<double[, ,]> R_matrix;
        private double[,] Lambda;
        private double epsilon; 

        public LR(int size, int numberdsc, double[, ,] A_disc)
        {
            this.epsilon = 0.01;
            this.system_size = size;
            this.discrete_num = numberdsc;
            this.Discret_matrix = new double[this.discrete_num, this.system_size, this.system_size];
            this.Lambda = new double[this.system_size, this.discrete_num];
            this.A_matrix = new List<double[, ,]>();
            this.L_matrix = new List<double[, ,]>();
            this.R_matrix = new List<double[, ,]>();

            for (int l = 0; l < this.discrete_num; l++)
                for (int i = 0; i < this.system_size; i++)
                    for (int j = 0; j < this.system_size; j++)
                    {
                        this.Discret_matrix[l, i, j] = A_disc[l, i, j];  
                    }
            this.A_matrix.Add(Listadder(A_disc));
        }

        public double[,] calculate_eigvalue()
        {
            int l =-1;
            do
            {
                l++;
                this.R_matrix.Add(Listadder());
                this.L_matrix.Add(Listadder());
                this.A_matrix.Add(Listadder());
                for (int k = 0; k < this.discrete_num; k++)
                {
                    Lmatrixinit(l, k);
                    RfirstLine(l, k);

                    try
                    {
                        Lfirstcolumn(l, k);
                    }
                    catch
                    {
                        return null;
                    }
                    Rsecondline(l, k);
                    try
                    {
                        Lsecondcolumn(l, k);
                    }
                    catch
                    {
                        return null;
                    }

                    RiLine(l, k);
                    try
                    {
                        LiColumn(l, k);
                    }
                    catch
                    {
                        return null;
                    }

                    AnextMatrix(l, k);
                    
                }
            }
            while (this.A_matrix[l+1][0, 0, 0] - this.A_matrix[l][0, 0, 0] > this.epsilon);
             
                 for (int i = 0; i < this.system_size; i++)
                     for (int k = 0; k < this.discrete_num; k++)
                         this.Lambda[i,k] = this.A_matrix[l][k,i,i];
            return this.Lambda;
        }

        private void AnextMatrix(int iter_num, int discrete_num)
        {        
            double currentsum = 0;

            for (int p = 0; p <= discrete_num; p++)
            for (int i = 0; i < this.system_size; i++)
                for (int j = 0; j < this.system_size; j++)
                {
                    currentsum = 0;
                    for (int n = 0; n < this.system_size; n++)
                        currentsum += this.R_matrix[iter_num][discrete_num-p, i, n] * this.L_matrix[iter_num][p, n, j];
                    this.A_matrix[iter_num + 1][discrete_num, i, j] += currentsum;
                }
        }

        private void LiColumn(int iter_num, int discrete_num)
        {
            double sum1 = 0;
            double sum2 = 0;
            double sum3 = 0;

            for (int ii = 2; ii < this.system_size; ii++)
            {
                for (int k = 1; k < this.system_size; k++)
                {
                    if (k < ii)
                    {
                        sum2 = 0;
                        sum1 = 0; 
                        for (int l = 1; l <= discrete_num; l++)
                        {
                            sum3 = 0;
                            for (int j = 1; j <= k - 1; j++)
                            for (int p = 0; p <= discrete_num; p++)
                                sum3 += this.L_matrix[iter_num][discrete_num - p, ii, j] * this.R_matrix[iter_num][p, j, k];
                            // this.L_matrix[iter_num][discrete_num, ii, 1] +=
                            sum2 += (this.A_matrix[iter_num][discrete_num, ii, k] - sum3) * this.R_matrix[iter_num][l, k, k];
                            //this.L_matrix[iter_num][discrete_num, ii, 1]
                        }
                        for (int j = 1; j <= k - 1; j++)    
                        for (int l = 0; l <= discrete_num; l++)
                        {
                            sum1 += this.L_matrix[iter_num][discrete_num - l, ii, j] * this.R_matrix[iter_num][l, j, k];
                        }
                        try
                        {
                            this.L_matrix[iter_num][discrete_num, ii, k] = (this.A_matrix[iter_num][discrete_num, ii, k] - sum1 - sum2) / this.R_matrix[iter_num][0, k, k];
                        }
                        catch (Exception ex)
                        {
                            throw ex;
                        }
                    }
                }
            }
        }

        private void RiLine(int iter_num, int discrete_num)
        {
            double sum =0;
            for (int ii = 2; ii < this.system_size; ii++)
            {

                for (int k = 1; k < this.system_size; k++)
                {
                    sum = 0;
                    if (k >= ii)
                    {
                        for (int j = 1; j <= ii - 1; j++)
                        { 
                            for (int l = 0; l <= discrete_num; l++)
                                sum += this.L_matrix[iter_num][discrete_num - l, ii, j] * this.R_matrix[iter_num][l, j, k];
                        }
                    }

                    this.R_matrix[iter_num][discrete_num, ii, k] = this.A_matrix[iter_num][discrete_num, ii, k] - sum;
                }
            }
        }

        private void Lsecondcolumn(int iter_num, int discrete_num)
        {
            double sum1 = 0;
            double sum2 = 0;
            double sum3 = 0;

            for (int ii = 2; ii < this.system_size; ii++)
            {
                sum2 = 0;
                sum1 = 0;
                for (int l = 1; l <= discrete_num; l++)
                {
                    sum3 = 0;
                    for (int p = 0; p <= discrete_num; p++)
                    sum3 += this.L_matrix[iter_num][discrete_num - p, ii, 0] * this.R_matrix[iter_num][p, 0, 1];
                    // this.L_matrix[iter_num][discrete_num, ii, 1] +=
                    sum2 += (this.A_matrix[iter_num][discrete_num, ii, 1] - sum3) * this.R_matrix[iter_num][l, 1, 1];
                   //this.L_matrix[iter_num][discrete_num, ii, 1]
                }

                for (int l = 0; l <= discrete_num; l++)
                {
                    sum1 += this.L_matrix[iter_num][discrete_num - l, ii, 0] * this.R_matrix[iter_num][l, 0, 1];
                }
                
                try
                {
                    this.L_matrix[iter_num][discrete_num, ii, 1] = (this.A_matrix[iter_num][discrete_num, ii, 1] - sum1 - sum2) / this.R_matrix[iter_num][0, 1, 1];
                }
                catch (Exception ex)
                {
                    throw ex;
                }
            }
        }

        private void Rsecondline(int iter_num, int discrete_num)
        {
            for (int ii = 1; ii < this.system_size; ii++)
            {
                for (int l = 0; l <= discrete_num; l++)
                    this.R_matrix[iter_num][discrete_num, 1, ii] += this.L_matrix[iter_num][discrete_num - l, 1, 0] * this.R_matrix[iter_num][l, 0, ii];
                
                    this.R_matrix[iter_num][discrete_num, 1, ii] = this.A_matrix[iter_num][discrete_num, 1, ii] - this.R_matrix[iter_num][discrete_num, 1, ii];
              
            }
        }

        private void Lfirstcolumn(int iter_num, int discrete_num)
        {
            for (int ii = 1; ii < this.system_size; ii++)
            {
                for (int l = 1; l <= discrete_num; l++)
                    this.L_matrix[iter_num][discrete_num, ii, 0] += this.A_matrix[iter_num][discrete_num - l, ii, 0] * this.R_matrix[iter_num][l, 0, 0];
                try
                {
                    this.L_matrix[iter_num][discrete_num, ii, 0] = (this.A_matrix[iter_num][discrete_num, ii, 0] - this.L_matrix[iter_num][discrete_num, ii, 0]) / this.R_matrix[iter_num][0, 0, 0];
                }
                catch (Exception ex)
                {
                    throw ex;
                }
            }

        }  

        private void Lmatrixinit(int iter_num, int discrete_num)
        {
            if (discrete_num == 0)
                for (int i = 0; i < this.system_size; i++)
                    this.L_matrix[iter_num][discrete_num, i, i] = 1;
        }

        private void RfirstLine(int iter_num, int discrete_num)
        {
            for (int i = 0; i < this.system_size; i++)
                this.R_matrix[iter_num][discrete_num, 0, i] = this.A_matrix[iter_num][discrete_num, 0, i];
 
        }

        private double[, ,] Listadder(double[, ,] source)
        {
            double[, ,] mat = new double[this.discrete_num, this.system_size, this.system_size];
            for (int l = 0; l < this.discrete_num; l++)
                for (int i = 0; i < this.system_size; i++)
                    for (int j = 0; j < this.system_size; j++)
                        mat[l, i, j] = source[l, i, j];
            return mat;
        }

        private double[, ,] Listadder()
        {
            double[, ,] mat = new double[this.discrete_num, this.system_size, this.system_size];
            return mat;
        }
    }
}
